Sampling the diffusion paths of a neutral vacancy in Silicon with SIEST-A-RT 
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We report a first-principles study of vacancy-induced self-diffusion in crystalline silicon. Starting 
form a fully relaxed configuration with a neutral vacancy, we proceed to search for local diffusion 
paths. The diffusion of the vacancy proceeds by hops to first nearest neighbor with an energy barrier 
of 0.40 eV in agreement with experimental results. Competing mechanisms are identified, like the 
reorientation, and the recombination of dangling bonds by Wooten- Winer- Weaire process. 

PACS numbers: 61.72.Ji, 71.15.Mb,71.15.Pd, 71.55.Cn, 
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I. INTRODUCTION 

Self-diffusion in homogeneous solids is a funda- 
mental process of mass transport, in addition, it is 
responsible of the annealing of implantation damage, 
crystal to amorphous transition and the nucleation 
of extended defects. The native defects, such as va- 
cancies and self-interstitial, have been identified as 
the prime entities controlling the self-diffusion, as 
well as the diffusion of impurities in semiconductor 
lattices. For example, in silicon, antimony impurities 
are vacancy diffusers and phosphorus are interstitial 
diffusers^. Therefore, self-diffusion in Si has been ex- 
tensively investigated by experimentsi^SiMi^ and simula- 

In covalently bonded silicon, the neutral vacancy is 
the most stable charge state near room temperature^*^. 
While the diffusion mechanism and the barrier for sili- 
con have been experimentally determined decades ago^, 
we still lack a detailed knowledge of the various acti- 
vated mechanisms associated with the neutral vacancy. 
With standard ab initio molecular dynamics methods, 
the number of events generated is always insufficient to 
give reliable data for migration energies. Accordingly, 
most of the focus has been on formation energies, distor- 
tion and reorientation geometries^i. Few papers reported 
a calculated activation energy and a migration path that 
are in agreement with the experimental dat a^^i'^^ . Until 
now, the only studies available for exploring the energy 
landscape were that of Munro and Wales^ who exam- 
ined the dynamics of a vacancy, among others, using a 
tight-binding approach, and the work of Kumeda et al. 
— using first-principles calculations. They identified the 
migration barrier and the underlying mechanism. How- 
ever, a complete study for the topology of the landscape 
around a stable minimum is still needed. 

The principal goal of this study is to give a com- 
plete description of the diffusion mechanism in elemen- 
tary semiconductors. As such, the SIEST-A-RT method 
is built up to develop an entire methodology for quantum 
mechanically precise and reasonably fast energy land- 



scape exploration method. Our simulations are per- 
formed on supercells containing 215 atoms. This super- 
cell size is computationally tractable and constitutes a 
reasonable size to simulate the host crystal. We generate 
the diffusion paths using the activation-relaxation tech- 
nique (ART)'"^-'^^, which can sample efficiently the energy 
landscape of complex systems. The forces and energies 
are evaluated using SIESTA2Si2S, a self-consistent den- 
sity functional method using standard norm-conserving 
pseudo-potentials and a flexible numerical linear combi- 
nation of atomic orbital basis set. Combining these two 
methods allows us to identify very efficiently diffusion 
paths of various energy height. 

This paper is constructed as follows. In the next sec- 
^^^^iWg^e discuss the SIEST-A-RT method; simulation de- 
tails are given in section Hill We apply the method to 
the vacancy in silicon, the corresponding activated dy- 
namics are detailed and discussed in section Hvl Finally 
we emphasize the capability of the method in identifying 
competing events in silicon and its possible extention for 
future applications. 



II. METHOD: SIEST-A-RT 

In order to sample the diffusion paths associated with a 
vacancy in bulk Si, we integrate the activation-relaxation 
technique (ART nouveau)^2iiS to the ab-initio program 
SIESTA. This approach is similar to the integration of the 
hybrid eigenvector-following method to a tight-binding 
and an ab-initio code as performed by Munro, Kumeda 
and Walesi^iSi as well as the dimer-method^ which was 
coupled with VASP^^ by Henkelman and Jonsson-^^i. 

As both ART nouveau and SIESTA are described else- 
where in the literature^SiiS, we review these methods only 
brieffy. 



A. ART nouveau 

The activation-relaxation technique is a method for 
sampling efficiently the energy landscape of activated sys- 
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tems.-2S^ Instead of following in details the thermal fluc- 
tuations, it searches directly for transition states; these 
are characterized as first-order saddle points in the en- 
ergy landscape. This method was refined recently to en- 
sure a controlled convergence to the saddle point (ART 
nouveau)4S An ART event can be decomposed into the 
following steps. Activation: (1) Starting from a local 
minimum, push the configuration in a direction chosen 
randomly in the 3N-dimensional space; stop when the 
lowest eigenvalue becomes negative or when the config- 
uration is far enough from the initial configuration. (2) 
Follow the eigendirection corresponding to this negative 
eigenvalue while minimizing the energy in the perpendic- 
ular directions; stop when the total force (parallel and 
perpendicular to the eigendirection) falls below a given 
threshold (0.1 eV/A). This is the transition state. Relax- 
ation: After pushing over the saddle point, minimize the 
energy using any standard minimization algorithm. 

In order to sample only the activated mechanisms as- 
sociated with the vacancy, we restrict the initial random 
projection for each event to the atoms surrounding the 
defect. Except for this initial selection, all atoms are 
allowed to move during the ART event. 

ART uses the Lanczos'''^ algorithm to extract the low- 
est eigenvalue and its corresponding eigenvector; 16 to 20 
force evaluations are sufficient, irrespective of the size of 
the system, to extract a stable eigenvalue value and eigen- 
vector. Recent unpublished tests indicate that the num- 
ber of force evaluations needed to reach a saddle point 
is about the same for the dimer method and ART nou- 
veau^, indicating no significant difference between the 
various activation-relaxation methods. 



B. Siesta 

Forces and energies are evaluated using SIEST A^^i^^ , 
a self-consistent density functional method using stan- 
dard norm-conserving pseudo-potentials of TrouUier- 
MartinSiii factorized in the Kleiman-Baylander form^. 
Matrix elements are evaluated on a 3D grid in real space. 
The one-particle problem is solved using linear combina- 
tion of pseudo- atomic orbitals (PAO) basis set of finite 
range. The quality of the basis can be improved by ra- 
dial fiexibilization, add more than one radial function 
within the same angular momentum (Multiple-^), and 
angular fiexibilization or polarization, add shells of dif- 
ferent atomic symmetry (different 1). 

Because ART requires a numerical derivative of the 
force to compute the curvature of the energy landscape, it 
is essential to obtain highly converged forces. The various 
parameters used are discussed below. 



C. Optimization of the parameters 

The introduction of a vacancy in the crystal affects 
considerably the relaxation of the simulation cell. Previ- 
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FIG. 1: "(Color on line)" Convergence tests for the energy 
per atoms as a function of Brillouin zone sampling for 64, 216 
and 512 cubic supercells. The lines are guide to the eye. (see 
the text). 



ous studies of the silicon vacancy using supercells rang- 
ing from 32 to 216 atoms with various Brillouin zone 
sampling techniquea^^i^^i^^i^^i^^i^^i^^i^'^'i^^ show a broad 
spread of formation energies (from 2.6 to 4.6 eV) and 
structural relaxation with {D2d, C2v, C'su, Td) symme- 
tries. This scattering in mainly due to various conver- 
gence problems like basis set size, Brillouin zone sam- 
pling and supercell size. A recent work of Probert and 
Payne^ propose a systematic methodology for accurate 
calculations of defect structure in supercells. We proceed 
in a similar way in using supercells with up to 512-atoms. 



1. Mesh cut-off and the basis set 

With the presence of a defect, the single-C basis (SZ), 
i.e., a single orbital per occupied state, is unsatisfactory: 
both the formation energy and reconstruction surround- 
ing the vacancy are far from experimental value. More- 
over, we find that the lack of overlap between the PAO 
induced discontinuities in the force, causing serious prob- 
lems to the application of ART. It is possible to address 
these problems by using the optimized basis generated 
using the procedure of Anglada et al^^ (SZ-optimized). 
Although much improved, the discontinuity in the force 
remains due to the presence of an underlying grid for 
the integration in real space which causes breaking of 
the translational symmetry. While the energy is con- 
verged with a real-space grid of about 0.03 A (a mesh 
cut-off of 35-40 Ry), we find that we need a grid space of 
at most 0.02 A (mesh cut-off of 50 Ry) to stabilize the 
force derivative; this spacing is comparable to the typical 
atomic displacement during the dynamics monitored by 
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ART. 



2. Brillouin zone sampling 

Because the defect is highly locahzed even in large su- 
percells, we need to have fully converged forces and to- 
tal energies for the system. In order to isolate Brillouin 
zone sampling effects, we study force and energy conver- 
gence with respect to the density of k-points. We use 
a Monkhorst and Pack''^ mesh to sample the Brillouin 
zone. Uniform meshes of shifted by (fco,fco,A:o) are 
used, with q going from 1 to 4 and ko = or 0.5. For 
q=l no offsets are used (fco — 0) to sample the T point. 
Offsets of fco — 0.5 are used for g = 2, 3, 4 for each cubic 
supercell size. 

Figure ^ displays the total energy per atom embed- 
ded in differently sized supercells with respect to q. The 
total energy is converged to ImeV/atom at a mesh of 
3^ for the 64 supercell, this corresponds to a density of 
0.031 A""'^ calculated using our LDA lattice parameter 
5.39 A. As discussed by Probert and Payne^^ this den- 
sity is sufhcient to converge the total energy for larger 
systems. It is equivalent to a mesh of 2^ for 216 and 512 
cells with densities of 0.031 A~^ and 0.023 A~^ respec- 
tively. Comparing the forces for unrelaxed vacancy using 
different meshes and supercells of 63, 215 and 511 atoms, 
we find that atomic forces converge faster than the total 
energy as a function of /c-point sampling: convergence 
to 1 me v/A per degree of freedom is reached with a 2'^ 
mesh for the 63 cell and with the P point for larger cells. 



defect-defect interaction, this component tends to vanish 
for the 215 supercell. 

Structural relaxation is also of a major importance, 
it informs us about the correctness of the atomic forces. 
The 63-atom supercell shows a wrong relaxation pattern. 
Tetrahedral (Td ) symmetry is a metastable minimum ir- 
respective of the density of k points used with a forma- 
tion energy of 3.99 eV ( see the second line on in tabled . 
Symmetry need to be broken artificially to get a Jahn- 
Teller distortion leading to a D2d symmetry with lower 
formation energy (3.59 eV). This is misleading and shows 
that size effects are important for supercells of 63 ionic 
sites and smaller. This size is insufficient to represent the 
system and the diffusion pattern correctly as discussed in 
more details in section llVDI 

For larger system, size effects tend to have minor effects 
on the relaxation geometry. Already at 215 cells, D2d 
symmetry is obtained even for low density of k-points (P 
only). The formation energy using P point for the 215 
supercell is underestimated by 6% compared to the best 
converged results we got with the 511 supercell, it tends 
to converge to the formation energy of larger models re- 
laxed with P-point only. 

We conclude, in agreement with previous worksiitSSiSi 
that a supercell of 215-atoms with P-point is the mini- 
mal size and sampling we can use in order to ensure that 
both energy and forces are sufficiently accurate to pro- 
vide results of a quality at least equal to experimental 
precision. 



III. DETAILS OF THE SIMULATION 



3. Supercell calculations 

As reported by Puska et alj^, supercell method has ob- 
vious drawback because of the interaction of the defect 
and its periodic replicas. If the defect-defect distance is 
not large enough, the electronic structure of an isolated 
defect is distorted because the deep levels in the band 
gap form energy bands with a finite dispersion. The size 
of the supercell restricts also the ionic relaxation. The re- 
laxation pattern is truncated midway between the defect 
and its nearest image. 

In order to identify the smallest cell-size acceptable 
for the simulation of diffusion mechanisms with negligi- 
ble size effects, we study the structural relaxation and 
the formation energies with respect to supercell size by 
employing supercells as large as 511 ionic sites. 

Formation energies are calculated for fully converged 
supercells with respect to the basis set and Brillouin zone 
sampling as detailed in sections III C II Hi C 21 Results for 
63 and 215 atoms supercells are overestimated compared 
to the 511-cell results by 9 % and 2% respectively as 
shown in tablej]. This error on vacancy formation energy 
for small supercell was previously reported by Puska et 
ali^ as due to the interaction of the vacancy with its 
periodic images and reflects a repulsive component of the 



We perform first-principles electronic structure calcu- 
lations based on the density functional theory within the 
local-density approximation (LDA)''*'"'^- for the exchange- 
correlation energy. The static zero-temperature calcula- 
tions for energies and forces are computed from SIESTA. 
Forces on the atoms are computed using the Hellman- 
Feynman theorem. The nuclear positions are optimized 
by using a conjugate gradient algorithm (CG) to mini- 
mize the total energy. Using the unit cell of two atoms 
and increasing the number of k points the equilibrium lat- 
tice constant of bulk Si converges to 5.39 A. The cell vol- 
ume and shape are kept fixed during the simulation. The 
crystalline 216 supercell is relaxed until the largest force 
component is of about 0.002 eV/A. One conjugate gra- 
dient (CG) step is necessary to relax the structure with 
vanishingly small force components. To create the start- 
ing configuration with a defect, one atom is removed from 
the ideal crystalline cell. The obtained structure contains 
a mono-vacancy characterized by a vacancy-vacancy dis- 
tance of 16.17 A, with tetrahedral symmetry (Td). The 
structure is then allowed to relax at constant volume. 
All the events are generated starting from a fully relaxed 
215-atom unit-cell with a single vacancy displaying an 
expected Jahn- Teller distortion: the four atoms around 
the vacancy are paired two by two with D2d symmetry. 



4 



TABLE I: Converged formation energies for (Ef) for different supercells and k-points meshes compared with various ab-miUo 
calculations. For each work we specify the DFT functionals used: local density approximation (LDA), screened-exchange 
LDA (sx-LDA) and conjugate gradient approximation (GGA) with the corresponding plane-wave (PW) energy cutoff when 
applicable. (TB-EF) denotes tight binding calculations with eigenvalue following method. In the D2d symmetry or the pairing 
mode, the distorted tetrahedron formed by atoms around the vacant site has two short bonds for atoms belonging to the same 
pair (atom-atom) distance and four long bonds between atoms from different pairs (pair-pair). The resulting symmetries of 
the defect are reported as well. Experimentally, Watkins et al^ measured a formation energy of 3.6 ± 0.5 eV and Dannefaer 
et al.^ 3.6±0.2 eV 



Authors 


Method 


Cell-size 


k-poirits 


Ef (eV) 


pair - 


Distances (A) 

pair atom - atom 


Symmetry 


This work 


LDA 


63 

215 

511 


33 
3=* 

r 

2^ 
2^ 


3.6 

3.95 

3.1 

3.36 

3.29 


3.34 
3.19 
3.34 
3.35 
3.35 


2.79 

3.19 

2.89 

2.9 

2.9 


Td 

D2d 
D2d 
D2d 


Lento and Nieminen— > 


LDA 

PW (15 Ry) 
sx-LDA 


32 

256 

256 


\ 4 ' 4 ' 4 

V 4 ' 4 ' 4 / 

V 4 ' 4 ' 4 


3.6 
3.515 


3.608 
3.490 


3.084 
2.950 


D2d 
D2d 

£'2d-outward 


Probert et alw~ 


GGA 

PW (12Ry) 


256 


2^ 


3.17 


3.72 


3.1 


D2d 


Puska et alw~ 


LDA 

PW (15 Ry) 


63 
215 


F 

2^ 
F 


2.86 
3.41 
3.27 


3.42 
3.38 


3.40 
2.89 


C2v 
D2d 
D2d 


Mercer et alw^ 


LDA 

PW (12-20 Ry) 


127 


2^ 


3.63 


3.50 


2.89 


D2d 


I^umeda et olw^ 


LDA 

PW(12 Ry) 
GGA 


63 

215 

215 


F 

F 
F 


3.32 

3.902 

3.336 








Seong and Lewis— 


LDA 

PW (10 Ry) 


63 
511 


F 

F 


3.65 
4.12 


3.53 


3.00 


D2d 


Munro and Wales — 


TB-EF 


63 
215 


F 
F 


3.73 
3.902 








Antonelli et alw~ 


LDA 

PW(16 Ry) 


63 


2^ 


3.49 


3.53 
3.52 


3.01 
3.40 


D2d 
D2d 



as described in section llVAl More than 120 events were 
started from this geometry, each launched in a different 
random direction, providing an extensive sampling of the 
energy landscape around the vacancy. The results from 
this search confirm a posteriori that this is indeed the 
minimum-energy configuration for the isolated vacancy. 
Using the set of parameters detailed in this section, one 
event vi^ith SIEST-A-RT requires about 700 force evalu- 



ations and takes about 3 days of CPU times on a single 
IBM Power 4 processor. 
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IV. RESULTS 
A. Neutral vacancy relcuxation 

The details of the relaxation around the vacancy de- 
pend on the local charge^°. Since there is one electron 
per dangling bond for the neutral charge state {V'^), it is 
expected to undergo a Jahn- Teller distortionSiS by form- 
ing two pairs out of the four dangling bonds of the va- 
cancy. This is what we see. During the relaxation of the 
215-atoni cell, the total energy of the system is reduced 
by about 1.71 eV, due to the distortion of the lattice and 
the relaxation of the atoms around the vacancy. This 
relaxation energy is larger than the 0.36 eV as obtained 
from a study on 63-atom cell by Seong et LewisiS and 
the 0.9 eV calculated on a HF cluster calculation^^ and 
compared to values from a more recent study— 1.186 eV 
using GGA on a 256 supercell. 

The relaxation can be described as follows. The va- 
cancy nearest neighbors (NN) first move away from each 
other by pair along the < 110 > axisi^. The net local 
displacement of 0.39 A per atom result into a tetragonal 
symmetry configuration {D2d)- The distorted tetrahe- 
dron formed by the atoms has four long bonds of 3.35 A 
and two short bonds of 2.9 A as shown in Fig.|3|^A). The 
formation of bonds between the two pairs weakens the 
back bonds which are then significantly stretchedi^, up 
to 2.49 A . The range of perturbation introduced by the 
defect in the 215 supercell affects up to 5 surrounding 
shells: the average amplitude of relaxation falls proges- 
sively from 0.39 A for the first shell to 0.074 A for the 
fifth shell which is of the order of atomic vibrations at 
room temperature. 

The amplitude and the type of displacement com- 
pare very well with the previous results using 
pprpii, 14^15,19,22,25,26^27^ rpj^g results appear to be 

at odds with recent screened-exchange LDA calculations 
on partially relaxed neutral Si vacancy, it indicates that 
the relaxation could be outward^^. More work remains 
to establish if LDA really fails here. Our results on the 
pairing mode are compared with some of these results in 
tabled 

The formation energy Ef is calculated by taking into 
account the distortion of the lattice, using the expression: 



Ef 



En-1 



N -I 



N 



-E 



N 



(1) 



Where N is the number of atoms and E]y is the re- 
laxed total energy of the ideal structure. Em-i is the 
total energy of the distorted system containing the de- 
fect. Within SIESTA, the formation energy of the mono- 
vacancy is 3.36 eV, which is in agreement with previous 
LDA and experimental studies. On the simulation side, 
Kelly et al.^ calculated 3.5 eV, Blochl et alM. 4.1 eV, 
Mercer et alM- 3.63 eV and Kumeda et al^ 3.32 eV. 
Experimentally, Watkins et alA measured 3.6 ± 0.5 eV 
and Dannefaer et al.^ 3.6±0.2 eV. A complete compari- 
son in reported in tabled 
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FIG. 2: " (Color online)" Top; (a) The total minimum-energy 
path for the simple diffusion. The saddle points identified in 
various ART events are indicated by the same labels (B,C and 
D) as in Table II. The path is generated by previous knowl- 
edge of the initial, the saddle point and final states, then the 
overall configuration is relaxed using the nudged-elastic-band 
method—. The contributions from diff'erent moving atoms 
are decomposed in two: (b) The path followed by the difi'us- 
ing atom; (c) Difltusion of the other atoms around the defect. 
Here The assymetry is due to reconstruction. 



A . 1 11 C 





FIG. 3: "(Color online)" The vacant site is visible in the 
snapshots A and C. A is the initial minimum and C is the 
ideal split vacancy site, the arrow shows the direction of the 
difi'usion. 



The energy of 511 supercell, Ej \ is found to be 
3.29 eV after a full convergence of the energy. Previous 
works by Wang et ali^ and Seong and LewisiS reported 
Ej^^= 4.12 eV and 4.10, respectively. A recent study 
by Probert et al. give GGA results of a fully converged 
calculation for a vacancy in BCC silicon lattice with 256 
atoms, they obtain a formation energy of 3.17 eV. 



6 



TABLE II: Interatomic distances along the minimum energy 
path for a simple diffusion. The positions correspond to the 
labels in Fig. I^Ja). A correcponds to the initial minimum 
with a relaxed vacancy: the active atom (1) is bonded to four 
neighbors, B,C and D are transition states. During the migra- 
tion bonds are stretched progressively until all the distances 
becomes equal at the configuration C. 



Distance (A) 


di-2 


di-3 


di-4 


di-5 


di-6 


di-T 


A 


2.49 


2.39 


2.39 


2.89 


3.34 


3.34 


B 


2.67 


2.67 


2.65 


2.65 


2.76 


2.73 


C 


2.68 


2.68 


2.68 


2.68 


2.68 


2.68 


D 


2.68 


2.87 


2.61 


2.69 


2.61 


2.69 




B. 


Simple 


diffusion paths 







Even a defect as simple as a neutral vacancy can be as- 
sociated with many activated mechanisms, some of which 
leading to migration: among the generated events, 60 % 
are simple diffusion mechanisms, while 15 % are reori- 
entation and the remaining are more complex events in- 
volving a pair or more of active atoms. 



1. Nearest neighbor diffusion 

The simplest migration process involves hopping of the 
vacancy to a nearest- neighbor site as shown in Fig. |2Ia) : 
a nearest-neighbor of the vacancy moves along the < 
111 > direction toward the vacant site (shown in Fig.|3\) 
and passes the metastable split-vacancy site before falling 
in the tetrahedral site previously formed by the vacancy. 
The energy barrier for this mechanism is 0.40 ± 0.02 eV. 
The error estimate comes from the convergence criterion 
to the saddle point as explained in section III Al This re- 
sult for the migration energy is in good agreement with 
the experimental findings of Watkins et ali^ who mea- 
sured 0.45±0.04 eV. In a work similar to this one, but 
using the CPMD code with GGA (BLYP) functional, 
Kumeda et alM-found a barrier of 0.58 eV. 

The diffusion path is reconstructed in detail using the 
nudged-elastic band methodS£ and relaxed until the to- 
tal force is of 0.2 eV/A. This path shows an unusually 
long transition plateau — ranging from dr=0.7 to 1.2 A — 
with a slight energy minimum at the split-vacancy site. 
At this point, the 3 back bonds are stretched and the dif- 
fusing atom forms weak bonds (length 2. 685 A) with the 
surrounding atoms of the six membered ring. From this 
site, a small displacement of 0.11 A is sufficient to over- 
come the barrier of 0.04 eV and complete the jump to the 
new minimum (see table . This value is of the order 
of the precision of the method; in addition, at room tem- 
perature this minimum would be smeared out by ther- 
mal vibrations. Because of a lack of a better approxima- 
tion, we use Vineyard's quasiharmonic approximations''^^ 



to get a rough estimate of the attempt frequency (see 
section HV H 3|) . approximating the transition plateau by 
a single barrier. 

Although the total diffusion path in Fig.|2Ia) is asym- 
metric, after separating the contribution coming from 
different moving atoms, we find that diffusion path is 
symmetric for the diffusing atom as shown in figure 
(Fig. Elb)). The asymmetry in the path shown in 
Fig. EIc) is due to the reconstruction: with four atoms 
reconstructing two-by-two, there are three possible ori- 
entations. Here, the orientation is different for the initial 
and final state. This should average out as the vacancy 
moves around. 

A comparison of the formation energies using the 
Stillinger-Weber^^ and the MEAM'^'' potentials for the 
vacancy and the split vacancy gives a migration energy 
of 0.43 and 0.37 eV, respectively. Although the formation 
energies for the vacancy are not well represented by em- 
pirical potentials the energy difference is in agreement 
with experiment and with our results. The minimum- 
energy path, however, is not well described by these po- 
tential and the local minima identified by Maroudas et 
al.'^^ during the migration are absent in our case. 



2. Reorientation process 

In addition to finding diffusion mechanisms associated 
with the vacancy, we also identify the trajectories re- 
sponsible for atomic rearrangement in the reconstructed 
configuration. As there are three pairing possibilities for 
the D2d symmetry, it is possible for the configuration to 
hop from one of these states to either of the two others. 

SIEST-A-RT generates easily this reorientation of the 
paired atoms. As illustrated in Fig. ^ starting with 
the reconstructed pairs formed initially — 1-4 and 2- 
3 — , atom "3" moves away form its neighbor along the 
< 111 > direction, breaks the weak bond with "2" and 
reaches the transition state, which is in the trigonal sym- 
metry (Csv). The relaxation continues until new bonds 
form between the pairs 1-2 and 3-4. The activation en- 
ergy for the reorientation of the neutral vacancy between 
two tetragonal orientations is 0.2 eV. Roberson et alm^ in 
a HF calculation scheme performed calculation on small 
clusters (44 silicon atoms) and found 0.37 eV; another 
cluster calculation, by Ogiit et alm^ found 0.32 eV. 

The deformation to the Csv can be observed experi- 
mentally. The experiment consists of applying a stress 
at high temperature in the dark, then cooling down to 
about 20K, relieving the stress, illuminating the sam- 
ple and monitoring the resulting vacancy species. The 
barrier measured from this alignment experiment is at 
0.23 eV, in excellent agreement with our calculations'^. 

Our simulations also show that the pairing mode can 
change during the migration of the vacancy. For the sim- 
ple diffusion plus reorientation, we find an overall energy 
barrier of 0.47 eV and a total displacement of 0.68 A at 
the saddle point. 
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Experiments under thermal equilibrium concentrations 
of native point defects have shown that the Si self- 
diffusion coefhcient D'g^ exhibits Arrhenius behavior, 



Doe 



(4) 



wiere Ha is a single activation enthalpy, Dq is the pref- 
axtov in a certain temperature range, T denotes the 
absolute temperature, and is the Boltzman's con- 
stant. Under the assumption that the diffusion proceeds 
rough discrete jumps of equal length r — the distance 
etween nearest neighbor equilibrium sites, in diamond 
cubic lattices r — where a is the lattice parameter 

the prefactor is written as24 



0,4 0.8 \2 

Total displacement 

FIG. 4: "(Color on line)" The reorientation path passing over 
the saddle point configuration. (I^S^F), the tetragonal axis 
is fixed throughout this process. At the saddle point the local 
atomic coordinates of the distorted lattice shows a trigonal 
symmetry. The pairing mode changes to one of the equivalent 
modes. 



Diffusion constants and diffusion rate 



Do 



6 



For" 



(5) 



w-lth To the jump, or attempt, frequency of the point 
defect from one equilibrium site to another, z denotes 
the number of neighbor sites, in the diamond lattice a 
vacancy can jump to one of the four neighbors. In the 
harmonic approximation, this attempt frequency is de- 
fined asSS 



3N 

n 



3N-1 



(6) 



The diffusion is a two step process, consisting of the 
creation and the subsequent migration of a vacancy. As- 
suming that both processes are thermally activated, the 
activation energy Ea for diffusion is the sum of the for- 
mation energy Ef and the migration energy Em- By 
collecting our data from the previous sections we get an 
activation energy of Ea = 3.5 eV, which is in good agree- 
ment with Ea = 4.14 eV recently measured by Bracht et 
alm^, and with the activation enthalpy for self-diffusion 



4.07 ± 0.2 eV obtained by Tang et 



usmg 



tight-binding MD. 

It is well known that there is a strong competition 
between the interstitial and the vacancy-mediated mech- 
anisms of self-diffusion. These two contributions were 
separated in a metal-diffusion experiments^ and the va- 
cancy contribution to the self-diffusion was measured to 
be 



(-4.14ey) 



-cm s 



Other experimental data-° shows that 



Cy'{T)d,{T) = O.Gexp 



{-4.03eV) 



-cm s 



(2) 



(3) 



Both data sets predict a self-diffusion dominated by va- 
cancies at low temperature, where Cy{T) is the concen- 
tration of vacancies at a certain temperature. 



where 



is) 

vl ' and v\ are the normal mode frequencies at 
the minimum and the saddle point respectively and the 
product does not include the imaginary frequency at the 
saddle point. The eigenvalues are computed by diagonal- 
izing the Hessian, obtained by numerical derivation with 
a step of 0.01 A. 

Accordingly, the migration entropy defined as the en- 
tropy difference between the saddle point and the equi- 
librium point configuration at low temperature can be 
written as 



AS- = 5: 



V 



/cs In 



3W-1 



,(«) 



(7) 



The phonon frequency corresponding the direction of the 
jump is removed from the numerator. In addition the self 
diffusion entropy evaluated experimentally is the sum of 
the diffusion entropy and the formation entropy Sy^ = 
Sy+Sy . Using the experimental Sy^ = 5.5kB by Bracht 
et alm^ and a 5/ = (5 ± 2)/cb as calculated by Blochl et 
al.^'^ and Dobson et al? , migration entropy is expected 
to lie within Ifcs and 2.5kB, our calculated data AS" = 
3.036fcB. 

The corresponding attempt frequency To = 8.65 x 
lO^^s"^ is similar to the result of Maroudas et ali^ 



(1.0539 X lO^^s"^), obtained with the StiUinger- Weber 
potential. The corresponding diffusion prefactor is Dq — 
3.141 X lO^^cm^s^^ . Our diffusion constant can be 
compared with results from Maroudas et al.'^^ by talc- 
ing into account the z factor that had been omitted, 
the resulting prefactor Dq ~ 3.88 x 10^^^cm^s^^ is in 
agreement with our results. Tang et al&^ report Dq = 
1.18 X 10~^cm^s~^ using TB-MD and thermodynamical 
integration, their migration energy (Em) is of 0.1 eV. 

Our numerical estimate of the diffusion rate T = 
Tqc^^/^'^'^ at room temperature gives a characteristic 
time for diffusion (r = l/F ) of the order of microsec- 
onds. This timescale is not directly accessible to ab-initio 
molecular dynamics, so simulations must generally be run 
at temperatures close to the melting point. With SIEST- 
A-RT, we can study the activated mechanisms from the 
local energy minimum, allowing us to identify accurately 
subtle mechanisms such as that associated with the split- 
vacancy site. 



C. Complex diffusion paths 

The SIEST-A-RT sampling of the energy landscape 
around the vacancy also identifies a number of high- 
energy barrier events which involve a complex atomistic 
rearrangements. Two types of such events occur. The 
first one can be described as a diffusion event mediated 
by bond exchange. In this case, the overall result is a 
simple vacancy hop but with a more complex transition 
state and a barrier of 2.4 eV. 

The second type of events involves what has been 
dubbed spectator eventsi^ (see Fig. where the di- 
amond network surrounding the vacancy is alterated. 
These spectator events all correspond basically to the 
bond-switching mechanisms proposed many years ago by 
Wooten, Winer and Weaire as mechanism for amorphiza- 
tion^. In recent years, the WWW mechanism has also 
been identified as the bond defect^^ and the four fold 
coordinated defect (FFCD)i&. 

Two types of spectator events are identified. In one of 
these events, two of the direct neighbors of the vacancy 
are pushed toward each other, leading to the saturation 
of the dangling bonds; as a consequence one of the long 
covalent bonds between the NN of the vacancy transform 
to a real bond in a form very similar to the bond defect 
complex identified by Marques et alii. The lowest barrier 
we found is 2.63 eV for a process involving a direct neigh- 
bor of the vacancy (one of the bonds is already weak), the 
final configuration is 2.16 eV higher than the initial min- 
imum. Kumeda et al. found the single spectator WWW 
barrier of 1.719 eV using LDA and 1.867 eV using BLYP 
functional 

The second type, atoms far from the vacant site are 
involved. The effect of the vacancy on the WWW event 
barrier is less strong for complex events located far from 
the vacancy without involving a direct neighbor. As the 
lattice distortion increases, so does the energy barrier. 
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FIG. 5: "(Color on line)" A spectator WWW process involving 
atoms far from the vacant site. Left panel: atoms in grey 
are the atoms involved in the exchange, the direction of the 
exchange is illustrated by the arrows. Right panel: the final 
configuration is a stable minimum by which a disorder in the 
crystalline lattice is introduced. 



For this mechanism, for example, the barrier ranges from 
3.31-3.8 eV, a value which can be compared to the I-V 
recombination process identified to be responsible for the 
amorphization of crystalline silicon of Marques et a!~ the 
barrier is estimated to be between 4.58 and 4.75 eV using 
ah initio calculation^. One can estimate the effect of the 
vacancy to lower the barrier by about 1.27 eV. Since the 
barrier is very high, this kind of events is unlikely to play 
any significant role in crystals at room temperature. 



D. Size effects 

As discussed, previously, it is necessary to use a large 
simulation cell in order to minimize the interaction be- 
tween the images of the vacancy. This is particularly 
noticeable in a 63-atom supercell, even though the defect- 
defect separation is larger than 4 inter-atomic distances. 
It is this interaction that cause the unrealistic relaxation 
pattern identified in previous studiesi^*2&, allowing the 
vacancy to relax both with a T2d and a D2d symmetries. 

Relaxing a 63-atom supercell with SIESTA and the pa- 
rameters described above, we find that two stable min- 
ima: the D2d (tetragonal) relaxation described previ- 
ously as well as a Td (tetrahedral) configuration, a sym- 
metry which is unstable in cells of 215 and 511 atoms. 
Within the later symmetry, all atoms relax inward, con- 
serving the tetrahedral symmetry of the system at equal 
distance of 3.4 A. The T2d structure can be reached with 
a single jump from the D2d site during a reorientation 
mechanism similar to that described in section IIV B 21 
crossing a small barrier of of 94 meV with an overall dis- 
placement of 0.64 A. The energy difference between the 
two symmetries is estimated to be 83 meV compared to 
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20 meV as obtained by Mercer et al. in a simulation 
also on a 63-atom celP^. We obtain qualitatively similar 
results with 2^ k-points mesh. 

The interaction between images also affects the bar- 
rier for the simple diffusion and the reorientation, which 
is underestimated by 50 % compared to the 215 cell re- 
sults. Interestingly, this underestimation is compensated 
by a higher formation energy, leading to very similar ac- 
tivation energies for the two size systems. 

The underestimation of the diffusion barrier in the 63- 
atom cell makes it easier for the vacancy to jump to the 
second-neighbor position. Although we saw many such 
events in the small cell, this mechanism was not gener- 
ated in the 215-atom box; there does not appear to be 
any direct path for second-neighbor diffusion in silicon. 

V. DISCUSSION AND CONCLUSIONS 

In this paper, we report on the study of the activated 
dynamics associated with a neutral vacancy in bulk crys- 
talline silicon. Following the work of Munro, Kumeda 
and Walesi^i^, we couple an ab initio code, SIESTA, 
with the activation-relaxation technique, ART nouveau, 
in order to sample the various activated mechanisms tak- 
ing place around the vacancy accurately and efficiently 
Defining moves directly in the energy landscape, it is as 
easy to generate the diffusion trajectory responsible for 
the high-energy bond-switching mechanisms than it is to 
find the path associated with the jump reorientation of 
the reconstructed state ( D2d ~* D2d)- 

Simulating more than 120 activated trajectories, we 
find a number of different diffusion mechanisms associ- 
ated with the vacancy. In particular, we recover the ba- 
sic vacancy hopping diffusion mechanism and show that 
the path associated with it includes a very metastable 
state at the split-vacancy site. Moreover, we can repro- 
duce without difficulty the jump between the symmetri- 



cally equivalent reconstructed states. We also generate 
a number of higher-energy diffusion mechanisms which 
reminiscent of those found in amorphous silicon as well 
as at the crystalline-amorphous interface. 

Because we obtain, with good accuracy, the transition 
state associated with each of these mechanisms, we also 
compute the transition rate in the harmonic approxima- 
tion. The result obtained is in good agreement with ex- 
periment as well as previous high-temperature TB-MD 
simulation's. 

Overall, various barrier and relaxation energies com- 
puted here are in agreement with results obtained previ- 
ously both experimentally and numerically. Using a large 
simulation cell and an extensive sampling of the energy 
landscape around the vacancy, we could also show that 
some of the previously obtained mechanisms were artifact 
due to the small size of the cell used. These results estab- 
lish the validity of the SIEST-A-RT approach as a unique 
tool for the study of the diffusion and relaxation mecha- 
nisms of defects in crystalline materials. This method can 
be easily extended to study semiconductor compounds, 
adatom diffusion, interlayer diffusion, all systems with a 
landscape too complicated for high-temperature molec- 
ular dynamics. Even in system as simple as a vacancy 
in c-Si, we found a number of activated mechanism as- 
sociated with non-trivial transition paths. The situation 
promises to be more complex in binary semiconductors, 
such as GaAs and GaN, where the bonding character and 
chemical composition are much more complex than for Si. 
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